
library(list)
library(doMC); registerDoMC()
library(random)
library(doRNG)

rm(list=ls())

iter <- 1

randSeed <- c(55762, 200542, 20002)

randSeed[iter] <- round(runif(1)*100000)

print(randSeed[iter])

set.seed(randSeed[iter])

load("data.RData")

treat.ISAFControl <- endorser[,1] == 1 | endorser[,1] == 2

dat <- as.data.frame(cbind(Y, indivcov, data[, c(colnames(villcov), colnames(distcov), "list.y", "list.treat")]))[treat.ISAFControl == TRUE, ]

names(dat)[1:4] <- paste("Q", 1:4, sep="")

dat$list.treat <- ifelse(dat$list.treat == "ISAF", 1, 0)

formula.temp <- as.formula(paste(" ~ ", paste(paste(colnames(indivcov), collapse = " + "), " + ", paste(colnames(villcov), collapse = " + "), " + ", paste(colnames(distcov), collapse = " + "))))

cov.mat.temp <- model.matrix(formula.temp, dat)

theta.start <- delta.start <- list()
lambda.start <- list()
x.start <- list()
for(i in 1:3) {
  theta.start[[i]] <- matrix(rnorm(4*29, mean = 0, sd = 5), nrow = 4, ncol = 29)
  theta.temp <- theta.start[[i]]
  lambda.start[[i]] <- matrix(rnorm(1*29, mean = 0, sd = 5), nrow = 1, ncol = 29)
  omega.start <- 1
  phi2.start <- 1
  delta.start[[i]] <- rnorm(29, mean = 0, sd = 5)
  x.start[[i]] <- rnorm(nrow(dat),
                        mean = as.numeric(cov.mat.temp %*% delta.start[[i]]),
                        sd = sqrt(9))
}

endorser <- endorser[treat.ISAFControl == TRUE,] - 1

MCMC <- 60000
burn <- 30000
thin <- 100

mle.estimates <- ictreg(as.formula( paste("list.y ~ ", paste(paste(colnames(indivcov), collapse = " + "), " + ", paste(colnames(villcov), collapse = " + "), " + ", paste(colnames(distcov), collapse = " + "))) ), J = 3, data = dat, method = "lm", treat = "list.treat", constrained = TRUE)

draws.list <- mvrnorm(n = 3, mu = coef(mle.estimates)[c(30:58, 1:29)],
                 Sigma = matrix(2, nrow = length(coef(mle.estimates)), ncol = length(coef(mle.estimates))))

procs <- 3

set.seed(45622)

res <- foreach(i=1:procs) %dorng% {

  endorse.options <- list(Y = list(Q1 = "Q1", Q2 = "Q2", Q3 = "Q3", Q4 = "Q4"),
                          data = dat, treat = endorser, covariates = TRUE,
                          formula = formula.temp,
                          delta.start = delta.start[[1]],
                          precision.delta = diag(1/9, 29), mu.delta = 0,
                          x.start = x.start[[1]], lambda.start = lambda.start[[1]], x.sd = FALSE,
                          tau.out = TRUE, MCMC = MCMC, burn = burn, thin = thin, seed.store = TRUE,
                          identical.lambda = TRUE, s.out = TRUE, mda = FALSE)
  
  result <- ictregBayes(paste("list.y ~ ", as.character(formula.temp)[2], collapse = ""),
                     data = dat, treat = "list.treat",
                     delta.start = draws.list[1, 1:29],
                     psi.start = draws.list[1, 30:58], burnin = burn, 
                     n.draws = MCMC, thin = thin, delta.tune = diag(.001, 29), psi.tune = diag(.0001, 29),
                     J = 3, constrained.single = "full", 
                     endorse.options = endorse.options, verbose = T, sensitive.model = "probit")
  
  return(result)
  
}

save(res, file = paste("results-combine.RData", sep = ""))

R.hat <- 2

while (R.hat > 1.1) { 

  burn <- 0

  res.last <- res
  
  res <- foreach(i=1:procs) %dorng% {
    
    N <- 1836
        
    endorse.options <- list(Y = list(Q1 = "Q1", Q2 = "Q2", Q3 = "Q3", Q4 = "Q4"),
                            data = dat, treat = endorser, covariates = TRUE,
                            formula = formula.temp,
                            precision.delta = diag(1/9, 29), mu.delta = 0,
                            x.sd = FALSE, update = TRUE, update.start = res.last[[i]],
                            tau.out = TRUE, MCMC = MCMC, burn = burn, thin = thin, seed.store = TRUE,
                            identical.lambda = TRUE, mda = FALSE, s.out = TRUE)
    
    result <- ictregBayes(paste("list.y ~ ", as.character(formula.temp)[2], collapse = ""),
                       data = dat, treat = "list.treat", burnin = burn, 
                       n.draws = MCMC, thin = thin, psi.tune = diag(.0001, 29),
                       J = 3, constrained.single = "full",
                       sensitive.model = "probit", endorse.options = endorse.options, verbose = T)
    
    result$x <- result$x[nrow(result$x), , drop = FALSE]
    result$s <- result$s[nrow(result$s), , drop = FALSE]
    result$beta <- rbind(res.last[[i]]$beta, result$beta)
    result$tau <- rbind(res.last[[i]]$tau, result$tau)
    result$lambda <- rbind(res.last[[i]]$lambda, result$lambda)
    result$omega2 <- rbind(res.last[[i]]$omega2, result$omega2)
    result$delta <- rbind(res.last[[i]]$delta, result$delta)
    result$psi <- rbind(res.last[[i]]$psi, result$psi)

    return(result)
        
  }
  
  save(res, file = paste("results-combine.RData", sep = ""))

}


